Brought to you by:

Kuiper Belt–like Hot and Cold Populations of Planetesimal Inclinations in the β Pictoris Belt Revealed by ALMA

, , , , , , and

Published 2019 March 7 © 2019. The American Astronomical Society. All rights reserved.
, , Citation L. Matrà et al 2019 AJ 157 135 DOI 10.3847/1538-3881/ab06c0

Download Article PDF
DownloadArticle ePub

You need an eReader or compatible software to experience the benefits of the ePub3 file format.

1538-3881/157/4/135

Abstract

The inclination distribution of the Kuiper Belt provides unique constraints on its origin and dynamical evolution, motivating vertically resolved observations of extrasolar planetesimal belts. We present ALMA observations of millimeter emission in the near edge-on planetesimal belt around β Pictoris, finding that the vertical distribution is significantly better described by the sum of two Gaussians compared to a single Gaussian. This indicates that, as for the Kuiper Belt, the inclination distribution of β Pic's belt is better described by the sum of dynamically hot and cold populations, rather than a single component. The hot and cold populations have rms inclinations of ${8.9}_{-0.5}^{+0.7}$ and ${1.1}_{-0.5}^{+0.5}$ degrees. We also report that an axisymmetric belt model provides a good fit to new and archival ALMA visibilities, and confirm that the midplane is misaligned with respect to β Pic b's orbital plane. However, we find no significant evidence for either the inner disk tilt observed in scattered light and CO emission or the southwest/northeast asymmetry previously reported for millimeter emission. Finally, we consider the origin of the belt's inclination distribution. Secular perturbations from β Pic b are unlikely to provide sufficient dynamical heating to explain the hot population throughout the belt's radial extent, and viscous stirring from large bodies within the belt alone cannot reproduce the two populations observed. This argues for an alternative or additional scenario, such as planetesimal being born with high inclinations, or the presence of a "β Pic c" planet, potentially migrating outward near the belt's inner edge.

Export citation and abstract BibTeX RIS

1. Introduction

The spatial distribution of dust in planetesimal belts around main-sequence stars provides significant constraints on the formation and dynamical evolution of planetary systems, both individually and as a population (e.g., Marino et al. 2018; Matrà et al. 2018a). Due to its vicinity and youth, the A6-type main-sequence star β Pictoris has long been considered a unique laboratory for exploring the outcome of the planet formation process. Over the course of more than 3 decades since the discovery of an infrared excess above the stellar emission (Aumann 1985), indicative of the presence of a dusty disk, numerous studies have been attempting to unravel the complexity of this planetary system.

At only 19.44 pc (van Leeuwen 2007), the disk's brightness and on-sky extent have made this system extremely suitable for resolved imaging. It was as early as 1984 when the belt's near edge-on geometry was discovered through optical coronagraphic observations (Smith & Terrile 1984), which showed starlight scattered off micron sized particles extending out to beyond 1000 au. Soon after, improved scattered light imaging of this system from the ground at optical and near-infrared (near-IR) wavelengths unveiled a complex disk structure with several asymmetries (Kalas & Jewitt 1995).

The near edge-on configuration makes β Pic's belt extremely well suited for studies of its vertical structure. Among the asymmetries, a warp in the spine of the on-sky brightness distribution within ∼80 au from the star was first unveiled through Hubble Space Telescope (HST) imaging and subsequently confirmed in all follow-up optical/near-IR scattered light observations (Mouillet et al. 1997; Heap et al. 2000; Golimowski et al. 2006; Boccaletti et al. 2009; Lagrange et al. 2012; Milli et al. 2014; Apai et al. 2015; Millar-Blanchaer et al. 2015). The presence of this warp was soon attributed to gravitational interaction with an unseen companion in the inner regions of the system (Mouillet et al. 1997; Augereau et al. 2001). It was not until 2009 that a giant planet, β Pictoris b, was discovered through direct imaging observations (Lagrange et al. 2009, 2010) at a projected distance of ∼9 au from the star (Lagrange et al. 2019), misaligned to the outer belt (Lagrange et al. 2012) in a way that is roughly consistent with the general expectation from the observed scattered light warp (e.g., Dawson et al. 2011; Nesvold & Kuchner 2015). Interestingly, the position angle of the CO and C i gas disks is similar to that of the inner disk seen in scattered light, suggesting a link between the two (Matrà et al. 2017; Cataldi et al. 2018), but no constraints have been reported so far on the presence of a warp in dust emission at longer wavelengths, due to a lack of angular resolution necessary to resolve this feature.

So far, the observed vertical distribution and width of the belt has received less attention. HST observations found that a sum of two Gaussians or Lorentzians best describes the observed vertical distribution of scattered light emission (e.g., Golimowski et al. 2006; Lagrange et al. 2012; Apai et al. 2015), with the scale height of the broad Gaussian component reaching values as high as ∼12 au at a distance of ∼80 au from the star. This corresponds to an aspect ratio of ∼0.15 at 80 au, similar to the aspect ratio inferred from modeling of polarimetric near-IR observations (∼0.137, Millar-Blanchaer et al. 2015).

The vertical distribution of planetesimal belts is a crucial observable as it reflects the inclination distribution of the constituent particles. In our own Kuiper Belt, the inclination distribution is broad and bimodal (e.g., Brown 2001; Kavelaars et al. 2009; Petit et al. 2011), and sets tight constraints on its dynamical history (Morbidelli et al. 2008), particularly on the outward migration of Neptune in the early stages of the Solar System's evolution (e.g., Nesvorný 2015). Therefore, accessing the inclination distribution of particles tracing planetesimals in extrasolar planetesimal belts is an attractive avenue to probe their dynamics and the effects of unseen planets.

With the assumption of an externally unperturbed and initially cold belt, the inclination distribution (particularly its width) could, for example, yield a direct link to the relative velocities of planetesimals. This allows us to probe the level of dynamical heating induced by unseen, large gravitational perturbers within the belt (e.g., Quillen et al. 2007). To date, only observations in optical/near-IR scattered light have had sufficient resolution to vertically resolve edge-on belts, such as around β Pic. Unfortunately, however, the radiation pressure felt by the small grains probed by these observations can significantly excite their eccentricities, leading to a consequential increase in their inclinations (Thébault 2009), which has prevented drawing a robust link between the observed vertical structure and the presence of large, unseen bodies required to create it.

In this paper, we tackle this issue using new ALMA interferometric observations that vertically resolve the β Pictoris belt at 1.33 mm. The observations are described in Section 2, and the main results are derived from image analysis in Section 3. In Section 4, we confirm these results by fitting the interferometric visibilities using parametric models to describe the belt's dust density distribution. In Section 5, we link the belt's observed vertical distribution to the inclination distribution of mm grains, and discuss its origin in the context of gravitational perturbations by large bodies within or exterior to the belt. Finally, in Section 6, we draw our conclusions and summarize our findings.

2. Observations

ALMA observations of the β Pictoris disk were obtained within its Cycle 2 (project code 2012.1.00142.S) using Band 6 receivers. Data were taken with the 12-m array in extended (one pointing centered on the star) and compact configuration (two pointings centered ±5'' from the star along the belt midplane), as well as with the Atacama Compact Array (ACA, one pointing centered on the star). Baselines covered a range between 9 and 1574 m, enabling us to be sensitive to structure between about 0farcs3 and 27''. A more complete description of the observations, including calibration, reduction, and imaging of CO and other lines, can be found in Matrà et al. (2017, 2018b). For continuum imaging, we flagged spectral channels around the frequency of the detected CO J = 2–1 line (230.538 GHz) and combined all the data sets from the different ALMA configurations, as well as the ACA. In doing so, we applied a constant rescaling factor to the visibility weights of each data set, where these factors (different for each data set) were determined from the scatter of the observed visibilities in Section 4.1. Then, we imaged the combined data set using the CLEAN algorithm (Högbom 1974) in multi-frequency synthesis mode. The data set covers ∼6.8 GHz of total bandwidth at an effective wavelength (frequency) of 1.33 mm (224.58 GHz). We used natural weighting of the visibilities to ensure maximum sensitivity while achieving a synthesized beam of size 0farcs32 × 0farcs27 and PA of −82fdg01; this implies a resolution of 6.1 × 5.3 au at the distance of the star (19.44 pc).

3. Results

3.1. New 1.33 mm Dust Continuum Data Set

Figure 1 shows the continuum emission map of the β Pictoris system obtained from the combined, naturally weighted, visibility data set. A peak signal-to-noise ratio (S/N) of 16 is reached at the location of the central star, for the rms noise level of 12 μJy/beam measured in an emission-free region of the continuum map. Emission from the belt is also clearly detected with a peak S/N of ∼14, as measured to the southwest (SW) of the star along the midplane. The spatially integrated 1.33 mm flux measured on the primary-beam-corrected map within a 16'' × 4'' box centered on the star is 20 ± 2 mJy (where the uncertainty is dominated by the 10% contribution from flux calibration; Fomalont et al. 2014).

Figure 1.

Figure 1. ALMA 1.33 mm CLEAN images of continuum emission from the β Pictoris belt obtained from the combined 12 m extended, compact, and ACA baselines, with natural weighting. The image is centered at the location of the star (detected here) and was rotated so that the belt midplane aligns with the horizontal axis, assuming the PA of the main disk observed in scattered light (29fdg3). The rms noise level achieved is 12 μJy beam−1 for a synthesized beam size of 0farcs32 × 0farcs27 (6.1 × 5.3 au, bottom left circle).

Standard image High-resolution image

The stellar contribution as measured from the peak flux at its location is 190 μJy, which should strictly be considered an upper limit due to extra disk emission lying along the line of sight to the star. In order to separate the stellar flux from the disk emission, we extract the subset of the data set containing visibilities at u − v distances >300. This u − v cutoff was chosen through inspection of CLEAN images, as an optimal trade-off to both ensure that all large-scale emission from the disk is filtered out, while retaining as much data as possible for maximum sensitivity. The obtained CLEAN image has an rms noise level of 18 μJy/beam for a beam size of 0farcs22 × 0farcs19, and a compact point source is detected at the stellar position. A simple 2D Gaussian fit within CASA yields a stellar flux of 80 ± 20 μJy at 1.33 mm. This detection is consistent with the expected stellar flux of ∼86 μJy obtained by extrapolating the Rayleigh–Jeans tail of the stellar blackbody-like emission from near-IR wavelengths. We therefore assume a stellar flux of 86 μJy at 1.33 mm throughout the rest of this paper.

3.2. Belt Radial Structure

In Figure 2, we show the profile of the continuum emission vertically integrated within ±20 au of the disk midplane, defined here as the position angle of the main disk observed in scattered light (29fdg3; Lagrange et al. 2012). Contrary to what is seen for CO emission in the same data set (green line), the millimeter-sized dust traced by the 1.33 mm continuum observations (blue line) shows no significant peak brightness asymmetry between the SW and northeast (NE) sides of the disk. This is in contrast with the SW/NE brightness asymmetry, between 30 and 80 au on either side the star, reported at 885 μm by Dent et al. (2014). The SW/NE ratio measured within the same radial distances as Dent et al. (2014), vertically integrated between ±20 au, is 1.03 ± 0.04 in the new 1.33 mm data set, consistent with no asymmetry. The same SW/NE ratio is 1.05 ± 0.03 in the 885 μm archival data set, indicating again a negligible (5 ± 3)% asymmetry in this radial region. As this measurement can be affected by imaging in the presence of incomplete u − v sampling and subsequent deconvolution via the CLEAN algorithm, we postpone further discussion of this issue in the context of visibility models to Section 4.6.

Figure 2.

Figure 2. Radial distribution of continuum 1.33 mm (blue), 885 μm (red), and CO J = 2–1 emission (green) obtained by spatially integrating images (such as Figure 1) within a height of ±20 au from the belt midplane. Each profile is normalized to its maximum to enable direct comparison.

Standard image High-resolution image

3.3. Belt Vertical Structure

Figure 3 (top) shows vertical profiles of emission perpendicular to the disk midplane as a function of distance to the star, which were obtained by integrating emission every 15 au along the radial (midplane) direction. As shown by comparing the profiles (blue shaded region) to the resolution element of our observation (black dotted line), the dust disk is clearly resolved in the vertical direction. Following the procedure described in detail in Section 3.2 and Appendix B of Matrà et al. (2017), we fit Gaussians (purple dashed lines) to the vertical profiles at each midplane location, and derive vertical best-fit Gaussian centroids (forming the disk spine, Figure 3 middle panel) and widths (measured as FWHM, Figure 3 bottom panel). We report two main findings as follows.

  • 1.  
    In contrast with the CO J = 2–1 line emission, we find that dust emission is not well represented by a single vertical Gaussian, being more centrally peaked and having broader wings extending out to ±40 au from the sky-projected midplane. This is evident in Figure 4, which shows the residual vertical profile obtained by integrating emission across the entire midplane after vertically centering it so that the centroid of the emission aligns with the midplane. The best-fit Gaussian profile (red, with a standard deviation of 8.2 ± 0.2 au) leaves significant residuals (Figure 4, bottom), whereas a double-Gaussian profile (green, with standard deviations of 5.1 ± 0.3 and 15.7 ± 1.4 au) considerably improves the fit, leaving no significant residual emission. We find that relative flux calibration systematics of ±10% between the data sets with compact and extended baselines (tracing vertically broader and narrower emission, respectively) cannot account for this discrepancy, which must therefore be physical in origin. This shows that the double-Gaussian or Lorentzian profiles needed to fit scattered light observations (Golimowski et al. 2006; Lagrange et al. 2012; Apai et al. 2015) reflect the vertical structure of the millimeter emission that traces the dust-producing planetesimals.
  • 2.  
    The belt spine (Figure 3, middle) presents little vertical displacement along the midplane and is largely consistent with the PA of the main disk observed in scattered light. This is significantly different from the CO J = 2–1 emission, whose spine presents an extra tilt angle dPA of ∼4° similar to that of the warp interior to ∼80 au seen in scattered light (Apai et al. 2015; Matrà et al. 2017). Therefore, we conclude that an inner disk tilt as large as observed in scattered light and CO observations is not detected in the structure of millimeter grains. Assuming the mm grains are tracing the planetesimal distribution, this shows that planetesimals are not perturbed into a warp—or at least not into one as pronounced as that observed in the distribution of the smallest grains observed in scattered light. At the same time, the tilt present in both the scattered light and the CO distribution, with the CO clump clearly offset by 4–6 au above the main disk midplane, suggests a connection between the two.

Figure 3.

Figure 3. Vertical structure of continuum 1.33 mm emission (blue) and CO J = 2–1 emission (red) from the belt. Shaded regions indicate ±1σ uncertainty ranges. Top: blue shaded regions represent normalized vertical emission profiles measured at midplane locations shown by the x-axis (see main text for details). Dotted lines represent the resolution element of the observations, whereas purple dashed lines are single-Gaussian fits to each of the observed profiles. Middle: shaded regions show the disk spines (i.e., the centroids of the vertical Gaussians fitted in the top panel) as a function of radius. CO J = 2–1 emission (red) is considerably misaligned with respect to 1.33 mm emission (blue) from the same ALMA data set. Bottom: vertical widths (FWHM) of the Gaussians fitted in the top panel, after deconvolving by the FWHM of the beam of our observations (projected on the perpendicular to the disk's midplane). These show similar heights for both CO and continuum emission.

Standard image High-resolution image
Figure 4.

Figure 4. Top: vertical profile integrated along the belt's sky-projected midplane of 1.3 mm emission (blue shaded region) as a function of height from the midplane. A best-fit Gaussian model is shown in red, whereas a best-fit double-Gaussian model is shown in green. Bottom: residuals obtained after subtracting the best-fit Gaussian (red) and double-Gaussian (green) models from the data, showing that the double-Gaussian model is significantly better at reproducing the observed vertical profile.

Standard image High-resolution image

4. Modeling

In this section, we aim to confirm our results inferred from images through radiative transfer modeling of the interferometric visibility data sets, enabling us to constrain the 3D distribution of mm dust in the belt while avoiding potential artifacts introduced by the imaging and CLEAN deconvolution processes.

4.1. Method

We model the β Pic belt using mm dust density distributions described in the following subsections, and the star as a point source with a flux density of 86 μJy located in the geometrical center of the circular belt. For each model realization, we solve the radiative transfer at 1.33 mm using RADMC-3D,8 assuming the dust radial temperature profile expected from a blackbody around a star of 8.7 L (Kennedy & Wyatt 2014). Since we are interested in the structure rather than the overall amount of material in the belt, we fix the dust opacity and fit for the belt mass as a normalization factor, noting that this does not affect our conclusions on the belt structure. In addition to the parameters that govern the belt dust density distribution and its overall mass, we fit for belt geometry parameters (inclination I and position angle PA).

Having created the model belt image, we produce visibility data separately for each of the ALMA arrays (12 m and ACA), for each configuration of the 12 m array and for each of the two mosaic pointings used for the compact configuration observations. We do so by first shifting the model belt image by R.A. and decl. offsets that we leave as free parameters in the fit, and multiplying it by the primary-beam response of the relevant array, at its appropriate sky location for each observation. We then Fourier transform the image to produce a grid of model complex visibilities, which we evaluate at the u − v locations sampled by the observation in question as described, for example, in Marino et al. (2016) and Tazzari et al. (2018). Finally, we add a point source representing the star to the model visibilities, located in the geometric center of the belt, and with its flux appropriately attenuated, given its location with respect to the primary-beam center.

The model visibilities are then fitted to the observed visibilities through a Monte Carlo Markov Chain (MCMC) method implemented through the emcee package (Foreman-Mackey et al. 2013). Eight nuisance parameters for the spatial shifts (in R.A. and decl., so two for each data set), two parameters to determine the viewing geometry (I and PA), one normalization parameter setting the mass of the belt for a fixed opacity, and n parameters governing the dust distribution make up the 11 + n-dimensional parameter space that we explore in our fit. Through the emcee package, we sample the 11 + n-dimensional posterior probability distribution of the explored parameters using the affine-invariant MCMC ensemble sampler of Goodman & Weare (2010). The posterior probability distribution is computed assuming uniform priors on all of the probed parameters, and a likelihood function proportional to ${e}^{-{\chi }^{2}/2}$, where χ2 is the usual definition of the chi-square function. The MCMC was run for each tested model with a number of walkers equal to 10× the number of free parameters and more than 2500 steps, ensuring each of the chains had reached convergence.

The uncertainty σ on each visibility data point is defined using the visibility weights w delivered by the ALMA observatory, where ${\sigma }_{{u}_{i}-{v}_{i}}=1/\sqrt{{w}_{{u}_{i}-{v}_{i}}}$. While assuming the delivered weights and hence the uncertainties are correct for different u − v points relative to each other, we left a constant rescaling factor, equal for all u − v points within a given data set, as a free parameter in our model fit of Section 4.2. This is justified by previous ALMA modeling studies also needing to rescale the weights to appropriately represent the true scatter in the visibilities (e.g., Marino et al. 2018). We find this rescaling factor to be different but very well constrained to 0.4508 ± 0.0006, 0.927 ± 0.014, 0.690 ± 0.003, and 0.698 ± 0.003 for the 12 m extended, ACA, and for the two pointings of the 12 m compact configuration, respectively. We therefore fix the rescaling factors to these values before producing the image in Figure 1 and do the same for all subsequent model fitting runs.

4.2. Radially and Vertically Gaussian, Axisymmetric Model with Radially Constant Aspect Ratio

We begin by assuming the simplest dust density distribution model with the minimum number of free parameters—namely a radially and vertically Gaussian, axisymmetric belt with radially constant aspect ratio. The distribution of the dust density ρ in the belt reads

Equation (1)

where r and z are cylindrical coordinates, rc and σ are the center and standard deviation of the radially Gaussian belt, h is the aspect ratio (here independent of radius), and ρ0 is a normalization factor proportional to the total dust mass in the belt, which we fit for.

The best-fit parameters of this model, defined as the ${50}_{-34}^{+34}$th percentiles of their marginalized posterior probability distributions, are listed in Table 1 (Model 1). Within the framework of this model, and in qualitative agreement with Figure 1, we find the disk to be near edge-on, with an inclination from face-on of ${86.6}_{-0.3}^{+0.4}$ degrees, and with its Gaussian radial distribution being centered at ∼105 au with a width of ∼92 au.

Table 1.  Best-fit Model Parameters

Parameter Unit Model 1 Model 2 Model 3 Model 4 Model 5
rc au ${104.9}_{-1.1}^{+1.1}$ ${107.5}_{-1.0}^{+1.0}$ ${106.3}_{-1.3}^{+1.4}$ ${107.8}_{-1.0}^{+1.0}$ ${105.2}_{-1.4}^{+1.4}$
σ au ${38.9}_{-1.1}^{+1.3}$ ${36.7}_{-1.3}^{+1.0}$ ${36.6}_{-1.2}^{+1.1}$ ${36.4}_{-1.3}^{+1.1}$ ${36.3}_{-1.2}^{+1.1}$
h   ${0.070}_{-0.004}^{+0.004}$ ${0.014}_{-0.006}^{+0.006}$ ${}^{* }{0.020}_{-0.008}^{+0.010}$ ${0.054}_{-0.004}^{+0.004}$ ${}^{* }{0.072}_{-0.012}^{+0.014}$
I ${86.6}_{-0.3}^{+0.4}$ ${88.6}_{-0.3}^{+0.2}$ ${88.8}_{-0.3}^{+0.4}$ ${89.3}_{-0.5}^{+0.4}$ ${89.5}_{-0.4}^{+0.3}$
PA ${30.0}_{-0.1}^{+0.1}$ ${29.7}_{-0.1}^{+0.1}$ ${29.7}_{-0.1}^{+0.1}$ ${29.8}_{-0.1}^{+0.1}$ ${29.7}_{-0.1}^{+0.1}$
h2   ${0.110}_{-0.006}^{+0.008}$ ${}^{* }{0.14}_{-0.02}^{+0.02}$
a   ${0.20}_{-0.04}^{+0.05}$ ${0.21}_{-0.04}^{+0.05}$
β   ${0.7}_{-0.2}^{+0.2}$ ${0.4}_{-0.2}^{+0.2}$
p   ${0.7}_{-0.2}^{+0.2}$ ${0.9}_{-0.1}^{+0.1}$ ${0.85}_{-0.07}^{+0.07}$
Free Parameters   14 16 17 15 16
${\chi }^{2}$   1456364.9 1456297.6 1456297.4 1456312.2 1456297.3
Δχ2   −67.3a −0.2b +14.6b −0.3b
ΔAIC   −63.3a +1.8b +12.6b −0.3b
ΔBIC   −38.9a +14.0b +0.33b −0.3b

Notes. Model 1: radially Gaussian ring with peak rc and standard deviation σ, with a single-component, Gaussian vertical profile with a radially constant aspect ratio h. Model 2: same as model 1, with a second, broader Gaussian component to the vertical profile, with aspect ratio h2. The ratio of the peak value of the second component to that of the first component is a. Model 3: same as model 2, but with a radially varying scale height for both Gaussian components, following H ∝ rβ, or equivalently h ∝ rβ−1. Model 4: same as model 1, but with a vertical profile described by a non-Gaussian functional form (where $\rho \propto \exp \left[-{\left(\tfrac{| z| }{\sqrt{2}{hr}}\right)}^{p}\right]$, see Equation (3)), with a radially constant aspect ratio h. Model 5: same as model 4, but with a radially varying aspect ratio h ∝ rβ−1. For Model 3 and 5, * symbols indicate that aspect ratios are quoted at a radius of 50 au.

aDifference with respect to Model 1. bDifference with respect to Model 2.

Download table as:  ASCIITypeset image

Figure 5, top left, shows the residuals obtained after subtracting the best-fit model visibilities from the data and imaging the residual visibility data set. The bottom left panel of Figure 5 shows the residual radial profile constructed as in Figure 2; the lack of significant residuals and SW/NE asymmetries reaffirms our conclusion that no clear departure from axisymmetry is present, within our measurement uncertainties, in the ALMA 1.33 mm data set.

Figure 5.

Figure 5. Top left: dirty image of the residual visibilities after subtracting the best-fit axisymmetric model with a vertical Gaussian density distribution (Model 1, see Section 4.2 and Table 1 for best-fit parameters) from the 1.33 mm data. Contours represent ±2, 4, 6 times the rms noise level. Bottom left: radial profile obtained from the residual image by vertically integrating within ±20 au of the midplane, as in Figure 2, showing no significant SW/NE asymmetry. Right: radially averaged vertical profile of the residual image, showing significant residuals analogous to the simple fit of Figure 4. These indicate that a single vertical Gaussian distribution is not a good fit to the data.

Standard image High-resolution image

On the other hand, the right panel shows the vertical profile obtained after radially integrating the belt's residual emission as done in Figure 4. This displays significant residuals very similar to those obtained by simply subtracting a single Gaussian from the vertical flux distribution in the image plane. This confirms that the vertical distribution of mm dust in the planetesimal belt around β Pictoris cannot be reproduced by a single Gaussian.

4.3. Double-Gaussian Vertical Structure Model, with Radially Constant Aspect Ratio

Motivated by the appearance of the vertical emission distribution in the image plane, we construct a belt model with a vertical distribution characterized by a narrower and broader Gaussian. The Gaussians are both centered in the belt's midplane (z = 0) and are characterized, respectively, by aspect ratios h and h2 leading to standard deviations (scale heights) hr and h2r at any given radius r. They comprise a fraction a and 1 − a of the total belt mass. The functional form of the dust density distribution thus reads

Equation (2)

This model is significantly better at describing the vertical distribution of the belt, as seen when comparing its residuals (Figure 6) to the single vertical Gaussian model (Figure 5). Comparing best-fit models (Model 2 and 1, respectively, in Table 1) through an Akaike or Bayesian Information Criterion (AIC, or BIC), which penalize models with a larger number of free parameters, clearly favors the model with the more complex vertical structure. Formally, the large BIC differential (or analogously the ratio of the model likelihoods) indicates the model with a double-Gaussian vertical structure is >108 times more likely to be better at describing the data than the model with a single Gaussian, providing strong evidence in support of our conclusion.

Figure 6.

Figure 6. Same as 5, but for a model where the vertical dust density distribution is the sum of two Gaussians (Model 2; see Section 4.3 and Table 1 for best-fit parameters). As demonstrated by the vertical profile of the residuals (right), this model produces a much better fit to the observations.

Standard image High-resolution image

The narrow Gaussian component has a vertical FWHM of ∼2 au at 150 au, which is similar to the size of our resolution element and only resolved through the forward modeling presented here. Note also that the inclination of the double-Gaussian model is significantly higher (closer to edge-on) than found for the single-Gaussian model. This is likely attributable to the single-Gaussian model being driven to a less edge-on configuration in an attempt to capture the broad sky-projected vertical distribution of the disk emission. The position angle east of north of the mm dust disk (${29.7}_{-0.1}^{+0.1}$) is marginally consistent within the uncertainties with the outer disk observed in scattered light (${29.3}_{-0.1}^{+0.2}$; Lagrange et al. 2012).

4.4. Adding Model Complexity: Radially Varying Aspect Ratio

The appearance of the residuals in Figure 6, after subtracting a model where the aspect ratio is the same at all belt radii, indicates that we have already achieved a good fit to the data. Nonetheless, we test whether we can detect a radial dependence of the vertical aspect ratio through a model where the scale height is a function of radius, and this dependence is parameterized through a flaring parameter β. The dust density distribution reads the same as in Equation (2), but with h and h2 being a function of radius following $h(r)={h}_{50\mathrm{au}}{\left(\tfrac{r}{50\mathrm{au}}\right)}^{\beta -1}$ and similarly ${h}_{2}(r)={h}_{\mathrm{2,50}\mathrm{au}}{\left(\tfrac{r}{50\mathrm{au}}\right)}^{\beta -1}$. Then, β = 1 indicates a constant aspect ratio h = H(r)/r as considered in the previous models, β > 1 a flared disk, and β = 0 a disk with a radially constant scale height.

The marginalized posterior probability distribution of β results in $\beta ={0.7}_{-0.2}^{+0.2}$, indicating that models with a scale height increasing with radius (β > 0) are significantly preferred, with a marginal improvement achieved for models with an aspect ratio decreasing with radius (β < 1). However, we find no evidence that the extra free parameter β leads to an overall improvement in the fit to the data, which is confirmed by visual inspection of the residuals. We therefore conclude that our data supports models where the scale height increases with radius, with a flaring parameter β consistent with one.

4.5. Functional Forms Other than Gaussians

Our results do not rule out that a different parametric prescription for the vertical density distribution could describe the data as well as the double-Gaussian parameterization adopted here. For example, we here attempt to use a model where the drop-off of the vertical density distribution differs from a Gaussian (as implemented by Ahmic et al. 2009). The functional form reads

Equation (3)

with $h(r)={h}_{50\mathrm{au}}{\left(\tfrac{r}{50\mathrm{au}}\right)}^{\beta -1}$ (as in Section 4.4), and $C\,={\left[{\int }_{-\infty }^{+\infty }{e}^{-{\left(\tfrac{| z| }{\sqrt{2}{hr}}\right)}^{p}}\right]}^{-1}$ ensuring normalization of the vertical distribution. Distributions with p < 2 allow the vertical distribution to drop faster than a Gaussian the further we move from the midplane, but to have broader wings; Ahmic et al. (2009) found this scenario to be preferred in scattered light observations of the β Pic belt, with p values around or below 1.

We tried two model runs with the above functional form and p left as a free parameter. In the first (Model 4), we fix the aspect ratio (β = 1) as done for the Model 1 and 2 runs. In the second run (Model 5), we also leave the flaring parameter β free. Our fitting results for both model runs, shown in Table 1, indicate that p ∼ 0.9, in broad agreement with the scattered light findings, and consistent with the broad-winged distribution derived from the image analysis.

We also find that best-fit Model 4 and 5 (with p = 0.9) are significantly better than a single Gaussian (Model 1). Compared with the double Gaussian, constant aspect ratio model (Model 2, see Table 1, bottom), we find that Model 4 (with β fixed to 1) is slightly worse, whereas Model 5 (with β ∼ 0.4) is just as good at reproducing the data, with extremely similar residuals. In summary, there is no evidence in the mm data that a model with a single population of mm grains that has a vertical drop-off rate that is shallower than Gaussian (Equation (3) with p ∼ 0.9, Models 4 and 5) is preferred to a double-Gaussian model (Model 2). In fact, the convergence of Models 4 and 5 to a broad-winged distribution confirms our main finding that there exists a population of high inclination (dynamically hot) particles in the β Pic disk. For later interpretation, we prefer to use the double Gaussian simply to facilitate direct comparison with the Kuiper Belt (see Section 5.2.3).

4.6. On the SW/NE Asymmetry in the mm Dust: Comparison with the 885 μm Data

Independent of the models explored previously, the vertically integrated radial profile of the residuals of the 1.33 mm data shows no significant SW/NE asymmetry (bottom left panel of Figures 5 and 6, where the shaded region is the ±1σ confidence interval). Averaging emission between 30 and 80 au on each side of the star, as in Section 3, we find a negligible residual flux difference between the NE and the SW of 82 ± 85 μJy. Varying the ranges over which we radially average does not affect this result. We conclude that we do not detect a significant SW/NE asymmetry in the new 1.33 mm data set.

Given the discrepancy with the 15% asymmetry reported in the 885 μm data set (Dent et al. 2014), we independently re-reduced the data set as described in Section 2.2 of Matrà et al. (2017). Just like the 1.33 mm, compact configuration data set, this consisted of a mosaic pointed ±5'' from the stellar location and along the belt's PA. We imaged the continuum using the tclean task (CASA v5.1.0) in multi-frequency synthesis, mosaic mode, using natural weights and multiscale deconvolution to best recover faint, extended emission. We obtained an image with synthesized beam size of 0farcs70 × 0farcs55 at a PA of 64fdg4, and measure an rms noise level of 70 μJy beam−1. The spatially integrated line flux measured on the primary-beam-corrected map as in the 1.33 mm data, is 66 ± 6 mJy.

As done for the 1.33 mm data, we extract the subset of the data set containing visibilities at u − v distances >300 to separate the stellar emission from the disk contribution. We find that the star is not significantly detected. This allows us to set a 3σ upper limit of <786 μJy on the stellar flux at 885 μm, which is consistent with the expected 196 μJy obtained through Rayleigh–Jeans extrapolation from near-IR wavelengths. We therefore include the star with a fixed flux of 196 μJy in the modeling that follows.

We simultaneously fit the visibilities obtained for both pointings as described in Section 4.1, using the single and double vertical Gaussian models described in Sections 4.2 and 4.3. For each of the two models, we fix the model parameters to the best-fit values found for the higher resolution and sensitivity, 1.33 mm data set (Table 1, Model 1 and 2). We leave as free parameters the scale factor ρ0 governing the total flux of the belt, which is different at the shorter wavelength, as well as an R.A. and decl. offset and a weight-rescaling factor for each of the two pointings that should be different from the 1.33 mm data sets.

Figures 7 and 8 show residual maps for the single and double vertical Gaussian models, respectively. We find the residuals, particularly the radial and vertical profiles, to be remarkably similar to those obtained for the 1.33 mm observations. In fact, despite the lower resolution and sensitivity, the 885 μm data set is also better fitted by a model with a double over a single-Gaussian vertical distribution, and also shows no significant SW–NE asymmetry. Averaging emission between 30 and 80 au on each side of the star, as in Dent et al. (2014), we find a negligible residual flux difference between the NE and the SW sides of 60 ± 130 μJy.

Figure 7.

Figure 7. Same as 5, but for a fit to the 885 μm data (see Section 4.6). Despite the lower angular resolution, the residuals show remarkable similarity to those obtained from 1.33 mm data.

Standard image High-resolution image
Figure 8.

Figure 8. Same as 6, but for a fit to the 885 μm data (see Section 4.6). Despite the lower angular resolution, the residuals show remarkable similarity to those obtained from 1.33 mm data.

Standard image High-resolution image

An interesting hypothesis is that the two Gaussian populations have different spectral indices, which could be caused by a difference in the size distribution of grains between the two populations. While this can in principle be tested with resolved observations at different wavelengths, such as presented here, the coarse resolution of the 885 μm data (which, at 100 au from the star, is of order the scale height of the broad Gaussian population) and the relatively small difference in wavelength between the two data sets do not allow us to test this hypothesis. Future ALMA observations at a S/N and resolution comparable to the 1.33 mm data set presented here, and at a wavelength as far as possible from 1.33 mm, are needed to explore this idea.

5. Discussion

In Sections 3 and 4 we modeled ALMA observations of the β Pictoris disk at 1.33 mm and 885 μm, which show the belt's vertical dust density distribution is much better fitted by the sum of two Gaussian distribution compared to a single-Gaussian distribution. In this section, we will link this vertical profile to the distribution of particle inclinations and connect it to their dynamical excitation to infer the potential action of large stirring bodies within or exterior to the belt.

5.1. Connecting a Belt's Vertical Density Distribution to the Planetesimals' Inclination Distribution

The observed inclination distributions of planetesimals in the Kuiper Belt was first investigated by Brown (2001, B01 hereafter). The authors used a Gaussian multiplied by a sine function to model the observed inclinations of Kuiper Belt Objects, and showed that this functional form naturally arises from N-body simulations from gravitational perturbations, or stirring, within a planetesimal disk with initially zero inclinations. The function reads $f(i){di}\propto \sin i\ {e}^{-\tfrac{{i}^{2}}{2{\sigma }_{i}^{2}}}{di}$ and for small angles corresponds to a Rayleigh distribution, $f(i)\propto \tfrac{i}{{\sigma }_{i}^{2}}{e}^{-\tfrac{{i}^{2}}{2{\sigma }_{i}^{2}}}$. This result is consistent with other dynamical studies also finding that gravitational perturbations of planetesimals in a thin disk lead to Rayleigh distributions of inclinations and eccentricities (e.g., Ida & Makino 1992; Lissauer & Stewart 1993). For direct comparison with the Kuiper Belt, we here follow the notation of B01 but note that dynamical studies typically express the inclination distribution as a function of the mean of the squares of the inclinations $\langle {i}^{2}\rangle $, with $\langle {i}^{2}\rangle =2{\sigma }_{i}^{2}$ in the expressions quoted previously. Figure 9 (left) shows an example inclination distribution (following B01) for a population of particles with a mean inclination of ∼7fdg9 (orange points) and shows how this is well matched to a Rayleigh distribution (blue line).

Figure 9.

Figure 9. Conversion from an inclination distribution of particles (left) to a measured latitude (center) and a scale height (right) above the belt midplane. Orange points show the inclination distribution used by Brown (2001) to reproduce the Kuiper Belt, how it is well matched (for the small angles considered here) to a Rayleigh distribution (blue line), and how it converts to a density distribution as a function of latitude and height above the midplane that is close to Gaussian (black dotted line).

Standard image High-resolution image

Under the assumption of circular orbits, this inclination distribution of planetesimal orbits f(i) in a belt at a given radius can be linked to the distribution of latitudes L(θ) (about a planetary system's ecliptic plane) at which the planetesimals are found (B01), through

Equation (4)

The integral can be simplified through use of the small angle approximation (sinx ∼ x for both i and θ), given the small inclination angles considered here. For example, Figure 9 (center and right) shows the latitudinal and height distribution derived from an inclination distribution (left) with a mean $\bar{i}$ of 7fdg9 and an rms of 8fdg9, corresponding to the hot population of the β Pic belt. At the height corresponding to 3σ of the distribution (16.5 au, implying a latitude of 18fdg9), the small angle approximation is accurate to within 2%. Taking the small angle approximation, the integral in Equation (4) becomes

Equation (5)

and can be solved analytically leading to the expression

Equation (6)

Given the small inclinations i expected to inevitably lead to small latitudes θ, we can again use the small angle approximation to obtain $\cos \theta \sim 1$ and $\mathrm{erf}\left(\sqrt{\tfrac{\tfrac{{\pi }^{2}}{4}-{\theta }^{2}}{2{\sigma }_{i}^{2}}}\right)\sim 1$ (the latter easily shown given, e.g., θ ≪ 5σi and small θ). Then, we come to the conclusion that the latitudinal distribution of particles with a Rayleigh distribution of inclinations (Figure 9, left) is simply a Gaussian with standard deviation equal to the σi parameter of the inclination distribution (i.e., $L(\theta )\propto {e}^{-\tfrac{{\theta }^{2}}{2{\sigma }_{i}^{2}}}$). We verify this in Figure 9 (center) by comparing an evaluation of the numerical integral in Equation (4) with the Gaussian derived analytically previously.

Finally, for small angles the distribution of heights (rather than latitudes) of the planetesimals will also be very well approximated by the Gaussian

Equation (7)

where ${\sigma }_{z}=H={hr}\sim {\sigma }_{i}r$. Thus, using a Gaussian vertical distribution in our model fits of Section 4 is justified by the expectation of a Rayleigh distribution of inclinations for a belt whose dynamics is dominated by gravitational stirring from bodies within the belt alone (e.g., Lissauer & Stewart 1993; Stewart & Ida 2000). The aspect ratio of an observed vertical Gaussian distribution can then be used to derive the mean $\bar{i}$ and root mean square ${\mathrm{rms}}_{i}=\sqrt{\langle {i}^{2}\rangle }$ of the inclination distribution of belt particles, following

Equation (8)

In Sections 3 and 4 we showed how the vertical distribution of mm dust in the β Pictoris belt is best described by a sum of two Gaussian distributions. The derivation provided then indicates that the corresponding distribution of particle inclinations is bimodal, corresponding to the sum of two Rayleigh distributions (Figure 10). From now on, we will refer to these as the hot and cold populations of particles, with high and low inclinations, respectively. The best-fit rms inclination derived from Model 2 is ${0.156}_{-0.009}^{+0.011}$ rad (${8.9}_{-0.5}^{+0.7}$ degrees) for the hot population and ${0.02}_{-0.01}^{+0.01}$ rad (${1.1}_{-0.5}^{+0.5}$ degrees) for the cold population.

Figure 10.

Figure 10. Same as Figure 9, but for a distribution of inclinations that is the sum of two Rayleigh distribution (left). The center and right panels confirm that this corresponds to a vertical density distribution that is the sum of two Gaussians (black dotted line), as observed for the β Pictoris belt in our data.

Standard image High-resolution image

5.2. On the Origin of the Inclination Distribution

5.2.1. External Stirring: β Pic b

The only known body that could contribute to structure in the belt is the directly imaged, super-Jupiter mass planet β Pic b. The most recent orbital determination for the planet indicates a best-fit semimajor axis of 8.8 au, an inclination to the line of sight Ip of 89fdg3 and a position angle Ωp of 32fdg1 (Lagrange et al. 2019). When comparing to the belt as resolved by ALMA, we confirm that the orbital plane of the planet and the belt are misaligned with respect to one another. This is thanks to our high-resolution measurement of the belt's position angle Ωb and inclination to the line of sight Ib from thermal emission, which for the first time are free of assumptions regarding the unknown phase function of the grains, as had affected previous determinations from scattered light observations.

The misalignment of the orbital plane of the planet and that of the disk can be characterized by the angle dI, the inclination of the planet orbit with respect to the disk midplane. A simple transformation between reference frames shows that

Equation (9)

We find that the belt's inclination Ib to the line of sight (see Table 1) for all the best-fit models to be consistent with that of β Pic b's orbit (Ip ∼ 89fdg3, but see Figure 5 in Lagrange et al. 2019, for an estimate of the uncertainty). If confirmed by future observations, this would indicate a planet-belt orbital plane misalignment dI ∼ 2fdg4, approximately equal to the difference in position angle between the belt and the planet's orbital plane.

Assuming this misalignment has been retained for the ∼23 Myr age of the system, the misaligned orbital plane of β Pic b should significantly affect the belt's vertical structure through secular perturbations affecting the belt's planetesimals. In particular, β Pic b should impose a forced inclination ip = dI ∼ 2fdg4 on the belt's particles and cause their inclinations to oscillate between 0 and 2ip (e.g., Dawson et al. 2011; Nesvold & Kuchner 2015). Particles further out in the belt are affected on longer timescales; this has long been considered the origin of the belt's warp observed in the inner regions out to ∼100 au from scattered light observations (Mouillet et al. 1997; Augereau et al. 2001), which trace the smallest grains in the collisional cascade.

In this scenario, we would expect (1) the warp to be also present for mm grains and (2) the vertical distribution of dust density, and hence mm emission, to be enhanced at 0 and 2ip where particles spend most of their time during their oscillations. An example of the expected disk morphology imparted by β Pic b is shown in Figure 1 of Nesvold & Kuchner (2015). This shows clear vertical flux enhancements displaced ∼ ± 9 au from the midplane at 100 au from the star, which are clearly ruled out by the disk spine and vertical profiles derived from the new ALMA data (Figure 3, central and top panel). Indeed, the mm continuum emission is centrally peaked and displaced by no more than 5 au vertically (at the 3σ level) throughout the disk midplane. This would limit any sky-projected warp, if present, to less than 2fdg8 (when measured at 100 au), which roughly corresponds to the forced inclination imposed by the planet.

In summary, we find that the disk spine and centrally peaked vertical emission profile are inconsistent with models of the secular perturbations imposed by β Pic b on the belt. Taken alone, this could indicate that (1) the planet is less inclined compared to the belt than currently believed; (2) another massive planet is present in the outer regions of the belt, dominating its dynamics; or (3) β Pic b has only recently been put on a misaligned orbit and has not had time to interact with the outer belt. Such an alternative picture would have to be reconciled with the scattered light and CO inner disk tilt, and with the hot and cold inclination populations reported here.

5.2.2. Gravitational Stirring within the Belt and Constraints on the Size of the Belt's Largest Bodies

The strongest constraint on the double-Gaussian vertical distribution of the belt comes from its radial outer edge, which is least affected by line-of-sight integration of emission originating at different orbital radii. Even if the belt is misaligned with β Pic b's orbit, this outer region has not had enough time to be affected by secular perturbations from the planet within the age of the system. Therefore, any structure in this region, and particularly the hot population, is either inherited from the gas-rich protoplanetary phase of evolution (e.g., Walmswell et al. 2013) or must be the result of dynamical excitation from a large body/bodies other than β Pic b.

We here consider gravitational stirring by large bodies within the belt itself. Encounters with such bodies will act as a source of dynamical heating, setting the velocity dispersion of particles, and producing a single, Rayleigh distribution of inclinations and eccentricities (see Kokubo & Ida 2012, and references therein). Given that the observed inclination distribution of the β Pic belt clearly deviates from a single Rayleigh distribution, we can rule out gravitational stirring from within the belt alone as the source of the observed vertical structure. Nonetheless, here we consider gravitational stirring as a potential source of either the hot or cold populations, while keeping in mind that other mechanisms are required to produce the other population.

Gravitational stirring leads to a relative velocity of particles of the form ${v}_{\mathrm{rel}}=\sqrt{1.25\langle {e}^{2}\rangle +\langle {i}^{2}\rangle }{v}_{\mathrm{Kep}}$, where vKep is the local Keplerian velocity (e.g., Lissauer & Stewart 1993). Since $\sqrt{\langle {e}^{2}\rangle }=2\sqrt{\langle {i}^{2}\rangle }$ (Ida & Makino 1992), the rms of the inclination distribution is a direct measure of the relative velocity of particles, following

Equation (10)

where h is unitless, r is in astronomical unit, M is in Solar masses, and the resulting vrel is in km s−1. Then, for an assumed stellar mass of 1.75 M (Crifo et al. 1997), the cold and hot populations in the β Pic belt (Model 2) indicate approximate relative velocities of 0.27 and 2.1 km s−1 at 50 au, dropping to 0.16 and 1.23 km s−1 at 150 au.

Viscous stirring from a large body (or bodies) will act to increase the velocity dispersion of planetesimals, and consequently their inclinations, over time. Following Ida & Makino (1993), assuming the system has been evolving from a perfectly cold disk (e, i ∼ 0 at t = 0 Myr) for the age of the system (t ∼ 23 Myr), we can then relate the measured inclinations to the product of the mass and the surface density of the largest bodies.

We assume we are in the dispersion-dominated regime and that stirring is dominated by large body/bodies rather than the small planetesimals themselves. Combining Equations (4.2) and (4.9) of Ida & Makino (1993) and solving for the inclination as a function of time, we obtain

Equation (11)

where the constant CI ∼ 2 (Ida & Makino 1993), Ω is the Keplerian frequency in rad s−1, M is the stellar mass in kg, t is time in s, and Σ is the surface mass density of the large stirring bodies (in kg m−2), each of mass M (in kg). The measured rms inclination of a planetesimal population then allows us to set a constraint on the mass times the surface density of the largest bodies of

Equation (12)

(where M is now in M, r is in astronomical unit, and t is in Myr). This leads to a value of 5.9 × 10−6 ${M}_{\oplus }^{2}$ au−2 for the cold population, and 0.022 ${M}_{\oplus }^{2}$ au−2 for the hot population, at 150 au in the β Pic belt. This corresponds to the purple lines in the Σ versus M space shown in Figure 11.

Figure 11.

Figure 11. Dynamical constraints on the mass and surface density of the largest bodies within the belt that could be producing (at 50 or 150 au) the cold or hot population of inclinations observed. The purple lines come from the observed hot and cold population inclination constraints (Equation (10)). The requirement that large bodies have not collided over the age of the system produces the surface density upper limits shown as the orange downward triangles. The intersection between the purple and orange lines represents the lower limit on the mass of the largest bodies. If the stirring bodies are massive enough that vesc ≫ vrel, this intersection will correspond to the requirement that their escape velocity equals the observed relative velocities of planetesimals in the belt (blue dashed–dotted line). Then, the presence of at least one large body within stirring distance of the radius in question produces the lower limit shown as the green upward triangles. The mass and surface density of the largest bodies are constrained to the portion of the purple line highlighted in green. A further constraint (not affecting the results here) is the lack of radial gaps in the belt's observed radial distribution. This implies that no gap of size equal to the chaotic zone has been opened over the lifetime of the system, setting a mass upper limit (brown leftward triangles in the upper panel). The red dotted line represents the surface density of the Minimum Mass Solar Nebula (MMSN).

Standard image High-resolution image

Producing the observed excitation however requires the bodies causing the stirring to not collide (and hence be destroyed) within the age of the system. The rate of mutual collisions for a population of bodies of equal mass M can be estimated as

Equation (13)

where σ = 4πR2 is the effective cross section of the collision (with R = (0.24M/ρ)1/3 being the radius of the large bodies, and ρ their material density), n is the number density of such bodies, and ${v}_{\mathrm{esc}}=\sqrt{{{GM}}_{\star }/R}$ is their escape velocity. As in Equation (8), we assume ${v}_{\mathrm{rel}}=\sqrt{6\langle {i}_{p}^{2}\rangle }$, where $\langle {i}_{p}^{2}\rangle $ is the rms inclination of the large bodies doing the stirring, which is likely primordial and here assumed to be 0.05. Setting the collision timescale (${t}_{{\rm{c}}}={R}_{\mathrm{col}}^{-1}$) to be greater than the age of the system, we can obtain an upper limit to the surface density of the large bodies given their mass, shown as the orange downward triangles in Figure 11. In log–log space, the slope of this surface density upper limit transitions from positive at the low mass end where vesc ≪ vrel, to negative at the high mass end where vesc ≫ vrel.

The intersection between the orange and purple lines indicates the minimum mass necessary for large bodies to produce the observed level of stirring but at the same time avoid potentially destructive collisions in 23 Myr. It can be shown analytically by combining Equations (9) and (11) at the high mass end (where vesc ≫ vrel) that requiring the collision timescale to be equal to the viscous stirring timescale leads to the planetesimals' relative velocities to be not exactly equal to but of order the escape velocity of the largest bodies (blue dashed–dotted line). This means that this minimum mass can be approximately estimated by simply setting the planetesimals' relative velocities to be equal to the stirring body's escape velocity.

A further requirement for the surface density of large bodies in the belt is that at least one large body is present at the radii where the stirred population is observed. In other words, if we consider the stirred population near, say, 150 au, there needs to be at least one large body such that its stirring zone encompasses 150 au. This sets a lower limit on the mass surface density (green upward triangles in Figure 11),

Equation (14)

where ${\rm{\Delta }}{r}_{\mathrm{stir}}=8\sqrt{3}r{[M/(3{M}_{\star })]}^{\tfrac{1}{3}}$ is the width of the stirring zone (equal to twice the large body's feeding zone, following Equation (3.4) of Ida & Makino 1993). The intersection of the purple line and green triangles represents the maximum mass for a body to stir a given population, and for at least one such body to exist so that its stirred region encompasses a given radius at which the population is observed.

Finally, we can consider the lack of observed radial gaps in the belt as a further constraint on the large body's mass. This is because a body sufficiently large will open a gap of width roughly equal to its chaotic zone (${\rm{\Delta }}{r}_{\mathrm{chaos}}\sim 3r{(M/{M}_{\star })}^{\tfrac{2}{7}}$; e.g., Wisdom 1980; Morrison & Malhotra 2015; Marino et al. 2018). Simply setting the width of the chaotic zone to be less than the resolution of our observations leads to an upper limit to the mass of the bodies within the belt. At 150 au, we find this mass to be <0.15 M. However, we need to consider the timescale necessary to open this gap. Following Equations (7) and (8) in Morrison & Malhotra (2015), the timescale for a 0.15 M body to open a gap at 150 au in the β Pic belt is ∼5.7 Gyr. Then, the lack of radial gaps can be caused not only by the chaotic zone not being resolvable, but also by the large body's gap-opening timescale being longer than the age of the system, where the latter sets a more conservative constraint on the maximum planet mass. At 150 au (50 au), we find that a stirring planet should be less massive than 75 (25) M to not produce observable gaps in 23 Myr (brown left-pointing triangles in Figure 11).

In reality, this estimate is uncertain because the clearing timescale of Morrison & Malhotra (2015) is defined as the "time required to reach 50% of the final survivor fraction within the cleared zone" and is therefore not to the timescale to produce an observable gap. The latter may be shorter, for example, if the observing sensitivity is sufficient to detect gap-belt surface brightness contrasts corresponding to a higher survival fraction of particles in the belt. At the same time, the edge-on viewing geometry would make it harder to detect such gap-belt contrast compared to a face-on belt, potentially making the timescale to produce a detectable gap longer. Detailed N-body simulations combined with radiative transfer and observing simulations are needed to accurately pinpoint this timescale (see, e.g., Marino et al. 2018), but are beyond the scope of this work.

Combining the dynamical constraints, we find that the mass of the largest bodies within the belt necessary to stir a hot population such as observed for β Pic over 23 Myr should be between 0.007 and 54 M at 150 au. On the other hand, if we consider viscous stirring as the origin of the cold population, we find that the largest bodies must have masses below 0.4 M at 150 au (or 0.15 M at 50 au) to have maintained such a low dynamical excitation. These masses are all far below and therefore consistent with current detection limits from direct imaging (Lagrange et al. 2018).

To conclude, we underline that viscous stirring by large bodies within the belt alone cannot produce the observed complexity of the vertical structure of the β Pictoris belt. However, it could produce either the hot or cold population, while requiring another mechanism to produce the other. If that were to be the case, the observed dynamical excitation of the hot (cold) populations could be produced through viscous stirring by bodies of mass 0.007–54 (<0.4) M at 150 au.

5.2.3. Outward Migration of a β Pic c

A bimodal distribution of inclinations such as observed in the β Pictoris belt was discovered for our own Kuiper Belt (KB) by B01 (after debiasing was taken into account; see also, e.g., Kavelaars et al. 2009; Petit et al. 2011). For both belts, we clarify that a single, broad non-Rayleigh distribution of inclinations cannot be excluded (e.g., Volk & Malhotra 2011); we just follow the original, bimodal functional form of B01 for direct comparison.

Figure 12 shows how the inclination distribution of the β Pic belt (namely, best-fit Model 2) compares with that of the KB from B01. Overall, we find lower inclinations for both the hot and cold inclination populations of the β Pic belt when compared to the KB, indicating an overall lower level of dynamical excitation. The fraction of all particles residing in the cold population are, however, consistent between the β Pic belt and the KB (${20}_{-4}^{+5}$% and ∼19%–26%), with both belts' masses being dominated by the hot population.

Figure 12.

Figure 12. Comparison between the possible inclination distributions as derived from our observations of β Pictoris (blue lines, Model 2) and the Kuiper Belt (KB, green and red lines).

Standard image High-resolution image

In the Kuiper Belt, the high inclinations seen for the hot classicals as well as the resonant and scattered populations set important constraints on its formation scenario (for a review, see, e.g., Nesvorný 2015). These populations likely formed in a massive planetesimal disk within 30 au and then got implanted into their present-day orbit by outward-migrating Neptune. The high inclinations can be attained through sweeping resonances (e.g., Levison & Morbidelli 2003) or scattering by Neptune during its migration (Gomes 2003). Nesvorný (2015) showed that the latter mechanism can work as long as Neptune's migration was slow (>10 Myr).

By analogy with the Kuiper Belt, we consider whether β Pic's structure could be produced by outward migration of a hypothetical β Pic c planet located at the inner edge of the belt. The presence of such a planet has been suggested to produce the SW/NE asymmetry observed for the CO gas (Dent et al. 2014; Matrà et al. 2017), as well as for the small dust grains at wavelengths below 24 μm (e.g., Telesco et al. 2005), following models predicting resonant planetesimals to concentrate at specific longitudes relative to the planet (Wyatt 2003). To produce the asymmetry observed for the CO, two conditions need to be met. First, the 2:1(u) resonance trapping probability (where the (u) indicates planetesimals in 2:1 resonance whose libration center increases from 180° during the migration) needs to be sufficiently high, and second, the planetesimal eccentricities need to be excited to a high enough level (we here require e > 0.3; see Figure 6, top row, of Wyatt 2003).

These eccentricities can be excited by resonant forces, causing them to increase as the radial extent of the migration increases (Equation (22) in Wyatt 2003). We find that the ratio between the semimajor axis of planetesimals after and before migration should be ≳1.8 for eccentricities to reach values above 0.3. For example, outward migration of a planet from 31.4 to at least 57.4 au (close to the inner edge of the belt as observed today) could have swept planetesimals in the 2:1(u) resonance from 50 to at least 91 au and produced sufficiently high eccentricities that an asymmetry may be observed. Given this migration must have taken place for at most the 23 Myr age of the system, we estimate a lower limit on the migration rate of 1.13 au Myr−1.

Ensuring the 2:1(u) trapping probability is sufficiently high (we here require it to be above 50%; see Equation (8) and Table 1 of Wyatt 2003) allows us to link this migration rate to the mass of the migrating planet. We find that a planet with mass of at least 28 M must have migrated at least ∼25 au to produce an asymmetry as observed for the CO. Therefore, it is possible that outward migration of a super-Neptune-sized planet may have produced the SW/NE asymmetry observed for the CO and for the small dust at short wavelengths. Can this, however, be reconciled with the lack of asymmetry for the mm-sized dust?

Wyatt (2006) showed that the appearance of resonant structure depends on the size of the observed grains, with grains below a critical size Dcrit falling out of resonance. This size depends solely on the migrating planet's mass (and some stellar parameters; see Equation (14) in Wyatt 2003). For a planet of at least 28 M, this size corresponds to 242 μm, although grains up to ∼10 times this size can show significantly smeared resonant structure. The observed emission in the ALMA 1.33 mm data set is dominated by grains of size similar to the observing wavelength or a few times smaller (see, e.g., Figure 6, bottom left in Wyatt 2003). It therefore seems plausible that the resonant structure could be smeared out to the extent that it is undetected in the 1.33 mm ALMA data as well as in the 24 μm (Telesco et al. 2005) and 70–500 μm (Vandenbussche et al. 2010) data sets. Observations at wavelengths even longer than presented here could clarify whether this is the origin of the lack of observed structure at 1.33 mm.

The model predicts the structure of the short-lived, smallest grains and the CO to be asymmetric. This means that the asymmetric resonance structure is only lost, or smeared out, at intermediate grain sizes (corresponding to mid to far-IR/mm wavelengths). Asymmetric resonant structure should therefore be seen for the smallest grains, for the CO, and for the large resonant grains of sizes ≫242 μm. Note, however, that while the distribution of CO and blow-out grains traces the collisional mass-loss rate from the belt (which is proportional to the relative velocity of particles multiplied by the square of their cross section per unit volume), the distribution of the large resonant grains simply traces the particles' cross sections per unit volume. This has the important implication that the SW–NE asymmetry will be much less pronounced and potentially hardly detectable for large resonant grains when compared to the blow-out grains and the CO. This means that the ALMA 1.33 mm and any longer wavelength data could probe resonant grains, but the resonant structure may be too faint for detection.

In addition to the SW/NE asymmetry and the high inclinations observed at the estimated location of the resonant clumps, the outward-migrating planet needs to be able to produce high inclinations all the way out to the outer edge of the disk (i.e., beyond the resonant clumps), as well as the tilt observed for the CO and blow-out grains. A hot population far out in the disk may be produced via outward scattering, as seen for the scattered population of the Kuiper Belt (Brown 2001). At the same time, the inner tilt seen for CO and small grains may be a projection effect caused by two asymmetric, azimuthal resonance clumps in front of and behind the plane of the sky in a non-perfectly edge-on configuration (see, e.g., Matrà et al. 2017, Figure 8, bottom).

Overall, this β Pic c migration scenario could qualitatively reproduce the observed features of the outer belt, but needs to be tested through future planet migration simulations. These should be tailored to reproduce β Pic's gas and dust structure as seen by ALMA, and should take into account the influence that β Pic b would have on such a planet, which we did not consider here.

6. Conclusions

In this work, we presented new ALMA observations of the β Pictoris belt at mm wavelengths. This provided us with the highest resolution picture of large, bound grains in the belt to date, allowing us to study their vertical distribution for the first time. We report the following findings, confirmed by both image analysis and modeling of the interferometric visibilities from the ALMA data.

  • 1.  
    The vertical distribution of mm grains cannot be fitted by a single Gaussian. In turn, this implies that a single Rayleigh distribution is not a good description of the particles' inclination distribution. We find a model with a bimodal distribution of inclinations, analogous to the Kuiper Belt, to be a much better fit to the data. rms inclinations are ${0.156}_{-0.009}^{+0.011}$ rad (${8.9}_{-0.5}^{+0.7}$ degrees) for the hot component, and ${0.02}_{-0.01}^{+0.01}$ rad (${1.1}_{-0.5}^{+0.5}$ degrees) for the cold component, with the hot component containing ${80}_{-5}^{+4}$% of the observed mass.
  • 2.  
    The position angle and inclination of the belt to the line of sight, here determined free of bias from the scattered light phase function, indicate significant misalignment between the belt's and β Pic b's orbital planes. Focusing on the belt's inclination to the line of sight, we find that it is consistent with that of β Pic b's orbit, which if confirmed would indicate that misalignment is maximum in or close to the plane of the sky.
  • 3.  
    After fitting an axisymmetric model to the new 1.3 mm and archival 885 μm visibilities, we find no significant evidence for the previously reported SW/NE brightness asymmetry in the continuum emission arising from mm grains. This is in agreement with previous 24–500 μm dust observations, but in stark contrast with the asymmetries seen in both gas and optical/near-infrared observations.
  • 4.  
    Similarly, we find no significant evidence for a tilt of the inner disk (i.e., a warp) as pronounced as observed in scattered light and CO observations. While the presence of a smaller tilt is not ruled out, this would now require a mechanism to produce a larger tilt in scattered light and CO gas with respect to the mm grains.

We then consider possible origins for the observed inclination distribution observed in the mm grains. We begin by considering the misalignment of β Pic b with respect to the belt, which should force the particles' inclination to oscillate between 0 and twice the planet's inclination with respect to the belt plane. This is ruled out by the lack of significant warping and by the vertical distribution observed by ALMA. Furthermore, the strongest constraint on the presence of the hot population comes from the belt's outer edge, which cannot be affected by β Pic b's secular perturbations within the age of the system.

We therefore explore alternative scenarios, such as viscous stirring from large bodies within the belt. The latter alone cannot explain the observed inclination distribution, but could be partly responsible for the excitation of either the cold or hot populations. We lay down a framework to constrain the mass and surface density of large bodies by connecting the observed inclination distribution to the expectation from viscous stirring. We find that producing the hot population would require large bodies of mass between 0.007 and ∼54 M at 150 au within the belt. Most of these bodies cannot produce observable gaps, as their gap-opening timescale is longer than the age of the system, although whether the most massive ones ≳10 M could is sensitive to the adopted "gap-opening" definition. On the other hand, keeping the low dynamical excitation of the cold population would require bodies to have masses below 0.4 M at 150 au.

Another potential explanation could come from the presence of an outward-migrating β Pic c at the inner edge of the belt. Such a planet could produce both the SW/NE asymmetry seen in CO and small grains through resonant sweeping, and explain the hot population of inclinations observed in a way that is analogous to the inferred outward migration of Neptune in the Solar System. A ≳28 M planet having migrated from ∼31 to 57 au in the system could produce the SW/NE asymmetry for the CO and small grains. This can likely be reconciled with the lack of a SW/NE asymmetry for grains observed between 24 μm and 1.33 mm, but requires tailored dynamical simulations (e.g., accounting for the influence of β Pic b on such a planet) to draw a more robust conclusion.

To conclude, new evidence from ALMA data combined with dynamical arguments indicates that β Pic b is unlikely to be the only large body responsible for the belt's observed vertical structure, unless planetesimals were born on high inclination orbits in the earlier protoplanetary disk phase of evolution. Instead, we have shown that the presence of a hot population of inclinations, analogous to the Kuiper Belt's, could be produced through dynamical excitation by additional large bodies interior and/or exterior to the belt. Detailed dynamical simulations combined with observations at longer wavelengths or higher resolution are necessary to confirm or rule out the scenarios proposed here.

L.M. acknowledges support from the Smithsonian Institution as a Submillimeter Array (SMA) Fellow. G.M.K. is supported by the Royal Society as a Royal Society University Research Fellow. This paper makes use of ALMA data ADS/JAO.ALMA#2012.1.00142.S. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada), NSC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ.

Facility: ALMA.

Footnotes

Please wait… references are loading.
10.3847/1538-3881/ab06c0